library(tidyverse) # for data visualization
library(metafor) # for meta-regressions
library(xtable) # for meta-regression table
library(broom) # for tidy table

#### 1. Load and clean data for meta-regressions #####

Data <- read.csv2("DataForMetaRegs.csv") # load data
Data <- escalc(measure = "RD", ai = n_vote_treatment, bi = n_no.vote_treatment, ci = n_vote_control,
               di = n_no.vote_control, var.names = c("effect_size", "variance"), data = Data)
Data$turnout_control <- Data$turnout_control / 100 # make turnout a proportion
Data$victory_marginpc <- Data$victory_margin / 100 # make margin of victory a proportion
Data$standard <- 0
Data$standard[Data$message_type == "standard reminder" |
                Data$message_type_2 == "standard reminder" |
                Data$message_type_3 == "standard reminder"] <- 1
Data$duty <- 0
Data$duty[Data$message_type == "civic duty" | Data$message_type_2 == "civic duty" |
            Data$message_type_3 == "civic duty"] <- 1
Data$past_turnout <- 0
Data$past_turnout[Data$message_type == "past turnout" | Data$message_type_2 == "past turnout" |
                    Data$message_type_3 == "past turnout"] <- 1
Data$neighbors <- 0
Data$neighbors[Data$message_type == "neighbors" | Data$message_type_2 == "neighbors" |
                 Data$message_type_3 == "neighbors"] <- 1
Data$advocacy <- 0
Data$advocacy[Data$message_type == "advocacy" | Data$message_type_2 == "advocacy" |
                Data$message_type_3 == "advocacy"] <- 1
Data$other <- 0
Data$other[Data$message_type == "other" | Data$message_type == "Hawthorne" |
             Data$message_type_2 == "other" | Data$message_type_2 == "Hawthorne" |
             Data$message_type_3 == "other" | Data$message_type_3 == "Hawthorne" |
             Data$message_type == "$10 incentive" | Data$message_type_2 == "$10 incentive" |
             Data$message_type_3 == "$10 incentive" | Data$message_type == "$5 incentive" |
             Data$message_type_2 == "$5 incentive" | Data$message_type_3 == "$5 incentive" |
             Data$message_type == "$2 incentive" | Data$message_type_2 == "$2 incentive" |
             Data$message_type_3 == "$2 incentive"] <- 1
Data$nonUS <- 1
Data$nonUS[Data$country == "United States"] <- 0
Data$squared_turnout_control <- Data$turnout_control ^ 2

#### 2. List descriptive statistics ####

table(Data$message_type)
table(Data$highest_office_level)
table(Data$primary)
table(Data$country)
table(Data$victory_margin <= 5)
table(Data$victory_margin > 5 & Data$victory_margin <= 10)
table(Data$victory_margin > 10 & Data$victory_margin <= 15)
table(Data$victory_margin > 15 & Data$victory_margin <= 20)
table(Data$victory_margin > 20 & Data$victory_margin <= 25)
table(Data$victory_margin > 25 & Data$victory_margin <= 30)
table(Data$victory_margin > 30 & Data$victory_margin <= 35)
table(Data$victory_margin > 35 & Data$victory_margin <= 40)
table(Data$victory_margin > 40)
table(is.na(Data$victory_margin))
table(Data$highest_office_election_system)
table(Data$highest_office)
table(Data$number_of_mails)
table(Data$number_of_countrywide_offices)

#### 3. Create 7 meta-regression models and put them in a table to test hypotheses ####

model0 <- rma(data = Data, effect_size, variance) # model with only intercept
model1 <- rma(data = Data, effect_size, variance, mods = ~ turnout_control) # model for H1
model2 <- rma(data = Data, effect_size, variance, mods = ~ turnout_control +
                squared_turnout_control)
model3 <- rma(data = Data, effect_size, variance, mods = ~ primary) # model for H2
model4 <- rma(data = Data, effect_size, variance,
              mods = ~ highest_office_level == "local") # model for H3
model5 <- rma(data = Data, effect_size, variance, mods = ~ victory_marginpc) # model for H4
model6 <- rma(data = Data, effect_size, variance, mods = ~ turnout_control + primary +
                (highest_office_level == "local")) # model for H1-H2-H3
model7 <- rma(data = Data, effect_size, variance, mods = ~ turnout_control + primary +
                (highest_office_level == "local") + duty + past_turnout + neighbors +
                advocacy + other + nonUS + number_of_countrywide_offices) # full model
model8 <- rma(data = Data, effect_size, variance, mods = ~ turnout_control +
                squared_turnout_control + primary + (highest_office_level == "local") +
                duty + past_turnout + neighbors + advocacy + other + nonUS +
                number_of_countrywide_offices)

convert_pval <- function(x){ # function to convert p-value into significance stars
  case_when(x < 0.001 ~ "***",
            (x >= 0.001 & x < 0.01) ~ "**",
            (x >= 0.01 & x < 0.05) ~ "*",
            (x >= 0.05 & x < 0.1) ~ ".",
            x >= 0.1 ~ "")
}

MetaRegDF <- full_join( # create tidy table
  broom::tidy(model0) %>% mutate(m = "model0") %>%
    mutate(term = "intercept") %>% mutate(my_text = paste0(
      round(estimate, 3), convert_pval(p.value), "\n(",
      round(std.error, 3), ")")) %>% select(term, my_text) %>%
    rename(`Model 0` = my_text), # model 0
  broom::tidy(model1) %>% mutate(m = "model1") %>%
    mutate(my_text = paste0(
      round(estimate, 3), convert_pval(p.value), "\n(",
      round(std.error, 3), ")")) %>% select(term, my_text) %>%
    rename(`Model 1` = my_text)) %>% # model 1
  full_join(
    broom::tidy(model2) %>% mutate(m = "model2") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 2` = my_text)) %>% # model 2
  full_join(
    broom::tidy(model3) %>% mutate(m = "model3") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 3` = my_text)) %>% # model 3
  full_join(
    broom::tidy(model4) %>% mutate(m = "model4") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 4` = my_text)) %>% # model 4
  full_join(
    broom::tidy(model5) %>% mutate(m = "model5") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 5` = my_text)) %>% # model 5
  full_join(
    broom::tidy(model6) %>% mutate(m = "model6") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 6` = my_text)) %>% # model 6
full_join(
  broom::tidy(model7) %>% mutate(m = "model7") %>%
    mutate(my_text = paste0(
      round(estimate, 3), convert_pval(p.value), "\n(",
      round(std.error, 3), ")")) %>% select(term, my_text) %>%
    rename(`Model 7` = my_text)) %>% # model 7
  full_join(
    broom::tidy(model8) %>% mutate(m = "model8") %>%
      mutate(my_text = paste0(
        round(estimate, 3), convert_pval(p.value), "\n(",
        round(std.error, 3), ")")) %>% select(term, my_text) %>%
      rename(`Model 8` = my_text)) # model 8

MetaRegDF$term <- c( # names of lines in meta-regression table
  "Intercept", "Turnout in control group", "Turnout squared",
  "Primary election", "Local election", "Margin of victory",
  "Message: civic duty", "Message: past turnout", "Message: neighbors", 
  "Message: advocacy", "Message: other", 
  "Not in USA", "Number of countrywide offices")
xtable(MetaRegDF) # table that can be exported to R Markdown
sink("PS_submission/table.html")
print(xtable(MetaRegDF), type = 'html') # print in HTML format
sink()

#### 4. Create plot for meta-regressions ####

PlotData <- escalc(measure = "RD", ai = n_vote_treatment,
                   bi = n_no.vote_treatment, ci = n_vote_control,
                   di = n_no.vote_control,
                   var.names = c("effect_size", "variance"),
                   data = Data)
MinMax <- data.frame(turnout_control = seq( # select range of predicted probabilities
  min(Data$turnout_control), max(Data$turnout_control),
  by = 0.01))
Pred <- data.frame(predict(
  model1, newmods = MinMax$turnout_control),
  turnout_control = MinMax$turnout_control)

Plot <- ggplot(Data, aes(x = turnout_control, y = effect_size)) +
  geom_line(data = Pred, aes(x = turnout_control, y = pred)) +
  geom_ribbon(data = Pred, aes(y = pred, ymin = ci.lb,
                               ymax = ci.ub), alpha = 0.3) +
  geom_point(aes(size = (0.5 / sqrt(variance)))) +
  scale_x_continuous("Turnout in control group",
                     limits = c(0, 1)) +
  scale_y_continuous("Effect of treatment",
                     breaks = seq(-0.05, 0.25, by = 0.05)) +
  scale_size(range = c(0, 5)) +
  geom_hline(aes(yintercept = 0)) +
  geom_vline(aes(xintercept = 0)) +
  theme(panel.grid.minor = element_blank(),
        strip.text.x = element_text(size = 15),
        strip.background = element_rect(fill = "white", colour = "black"),
        axis.title = element_text(size = 13.5),
        axis.text = element_text(size = 13.5),
        panel.background = element_rect(fill = "white"),
        legend.position = "none")
ggsave("Figure2.png", dpi = 1200) # save plot
View(ggplot_build(plot = Plot)$data[[2]]) 
